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14.  ABSTRACT 


This  report  results  from  a  contract  tasking  Middle  East  Technical  University  as  follows:  The  aim  of  this  work  is  to  gain  an  understanding  of  the 
unsteady  aerodynamics  of  flapping  motion  in  order  to  implement  the  data  obtained  to  the  design  of  a  Micro  Air  Vehicle.  The  second  task  is 
increase  the  experimental  capability  of  the  existing  2D  flapping  motion  experimental  setup  to  3D. 

The  main  objective  of  the  previous  works  was  to  understand  the  physics  and  the  aerodynamics  of  lift  generation  by  flapping  motion  during 
hover.  The  initial  phase  of  the  study  was  limited  to  2D  analysis  where  a  symmetrical  NACA  0012  airfoil  is  moved  upstroke  and  downstroke 
with  variable  speeds.  This  research  will  be  a  continuation  of  the  research  work  undertaken  by  D.  Funda  Kurtulu0,  during  her  Ph.D  studies  in 
ENSMA  Poitiers  France  since  2002.  With  this  project  it  is  intended  to  improve  and  enhance  the  experimental  studies  in  3D. 

Two  different  methods  of  analysis  are  used  during  this  study. 

1 .  Numerical  methods  (DNS  and  aeroelastic  simulations) 

2.  Experimental  techniques 

Expected  results:  The  results  of  the  project  will  be  a  part  of  the  NATO  AVT-149  Technical  Team  group  and  will  be  published  in  the  final  report. 
Numerical  Method: 

In  numerical  analysis,  2D  flapping  airfoil  and  3D  flat  plate  rigid  platforms  will  modelled  by  Direct  Numerical  Simulation  technique  (DNS)  using  a 
computer  programme,  based  on  an  available  commercial  code  Star-CD  (or  any  other  3D  unsteady  codes  that  will  be  decide  in  NATO  AVT  149 
technical  group)  for  low  and  high  Re  numbers.  For  the  first  trials,  3D  results  will  be  limited  to  1  degree  of  freedom  (rotation).  The  major  aim  of 
flapping  motion  research  is  based  on  the  understanding  of  the  relation  between  the  temporal  and  the  spatial  changes  of  the  wake  structure 
and  the  resulting  instantaneous  aerodynamic  forces  over  the  flapping  wings.  The  essential  physics  of  non-steady  airfoil  problems  can  be 
observed  from  simplified  two-dimensional  experiments,  and  the  interpretations  of  the  behaviour  can  be  supported  by  theoretical  or  numerical 
models.  In  addition  to  the  instantaneous  aerodynamic  forces,  pressure  distributions  and  vorticity  contours,  the  average  lift  and  drag  coefficient 
values  are  also  calculated.  The  3D  platforms  (flat  plate)  will  be  also  investigated  with  aeroelastic  simulations  using  MSC  Nastran  program 
which  uses  Doublet  Lattice  Method  for  aerodynamic  analysis  and  couples  with  its  structural  model.  The  program  could  handle  also  gust 
responses.  Some  preliminary  tests  are  performed  with  static  analysis. 


Experimental  Techniques: 

The  experimental  technique  is  based  on  the  flow  visualisation  studies  of  a  two-dimensional  wing  in  a  water  tank,  whose  motion  is  controlled 
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Abstract 

The  numerical  and  experimental  investigations  on  the  characteristics  of  low  Re  number 
regime  are  increasing  due  to  the  advances  in  micro-technologies  enabling  the  development  of 
Micro  Air  Vehicles  (MAV’s).  One  of  the  main  objectives  of  MAV  applications,  i.e.  constant 
position  surveillance,  reveals  the  need  to  focus  the  researches  on  hover  mode.  There  are  three 
generations  of  MAV’s  namely  fixed  wings,  rotating  wings  (like  helicopters)  and  wings  based 
on  micro  technologies  (MEMS,  flapping  or  vibrating  wings).  The  definition  employed  in 
Defense  Advanced  Research  Projects  Agency  (DARPA)  program  limits  these  crafts  to  a  size 
less  than  15  cm  in  length,  width  or  height.  This  physical  size  puts  this  class  of  vehicle  at  least 
an  order  of  magnitude  smaller  than  any  UAV  developed  to  date. 

The  studies  on  flapping  motion  flight  can  be  classified  into  two  main  parts  as  the 
zoological  configurations  and  the  simplified  configurations.  The  zoological  configurations  are 
performed  based  on  the  study  of  insects  or  birds.  Comprehensive  reviews  of  the  biological 
flight  mechanisms  could  be  found  in  Nachtigall  (1974);  Rayner  (1979,  1985);  Ellington 
(1984);  Norberg  (1985);  Azuma  et  al.  (1985);  Pennycuick  (1988)  and  Dudley  (1998). 

The  simplified  configurations  are  mostly  the  works  based  on  the  aerodynamics  of  the 
flow.  The  models  are  simplified  such  that  the  real  insect/bird  wing  geometries  are  replaced 
with  different  pre-defined  aerodynamic  profiles  (Shyy  et  al.  1999;  Ramamurti  and  Sandberg 
2002;  Platzer  and  Jones  2006).  Although  numerical  methods  are  widely  used  for  unsteady 
aerodynamic  problems,  it  is  highly  difficult  to  solve  the  full  3D  Navier-Stokes  equations  for 
complex  flows  like  the  ones  of  the  flapping  insect  wings.  Hamdani  and  Sun  (2000)  simulated 
a  series  of  impulsive  starts  at  different  accelerations  around  a  2D  insect  wing.  The  mean 
streamwise  velocity  field  of  the  wake  of  a  NACA  0012  airfoil  oscillating  in  plunge  at  zero 
ffeestream  velocity  and  at  a  zero  angle  of  incidence  at  the  neutral  position  was  calculated  by 
Lai  and  Platzer  (2001).  The  vortical  flow  patterns  in  the  wake  of  a  NACA  0012  airfoil 
pitching  at  small  amplitudes  are  studied  by  Koochesfahani  (1989)  in  a  low-speed  water 


channel  by  considering  the  effect  of  both  sinusoidal  and  non-sinusoidal  shape  of  the 
waveform.  The  aeroelastic  simulations  of  the  problem  are  another  important  concern  of  the 
researchers  in  this  domain.  Some  first  simulations  with  doublet  lattice  method  will  be 
represented  for  this  project. 

The  objective  of  the  present  work  could  be  classified  as  follows: 

1 .  Numerical  simulations 

a.  3D  Rigid  Platform  (Flat  plate)  1DOF  with  low  Re  number  simulations  (translational 
impulsive  start,  rotational  impulsive  start,  parametrical  study) 

b.  2D  Rigid  Platform  2DOF  high  Re  number  simulations  with  SD  7003  airfoil 

2.  Experimental  results 

a.  3D  Platform  translational  impulsive  start,  rotational  impulsive  start  visualization  in 
water  tank  with  a  parametrical  study 

b.  3D  flapping  motion  in  low  Re  number  (Re<2000)  by  using  existing  experimental  setup 
and  investigating  the  3D  effects  by  removing  one  of  the  end  plates. 

3.  The  aeroelastic  numerical  calculations 

a.  3D  flat  plate  methods  could  be  calculated  with  MSC.  Nastran  aeroelastic  module. 
Conventional  MSC.Nastran  Structural  Model  performs  Static  Analysis  and  Normal  Modes 
Analysis.  MSC.Nastran  Aerodynamic  Analysis,  uses  flat  plate  methods  (Doublet  Lattice 
Method).  The  program  maps  Aerodynamic  Forces  to  Structure  and  Maps  Structural 
Deflections  to  the  Aerodynamics 

b.  Unsteady  Gust  analysis  could  be  performed  with  this  module  to  the  3D  flat  plate 


1.  Computational  Method 


The  solutions  have  been  performed  in  two  aspects.  One  is  the  low  Reynolds  number  laminar 
solutions  of  flapping  profiles  in  hover  mode  and  the  second  part  is  related  to  the  results  of 
AVT  149  Technical  Team  pure  plunge  and  pitch/plunge  case  studies  around  SD  7003  profile. 
Additionally  different  sinusoidal  motions  at  low  Reynolds  number  has  also  been  studied  and 
compared  with  the  experimental  data  in  literature. 


The  low  Reynolds  number  results  for  hover  mode  were  given  in  Report  1  (Appendix  A)  and 
they  are  not  repeated  in  this  second  report.  The  additional  works  carried  out  from  January 
2009  are  summarized  in  this  Report  2. 

1.1  AVT  149  Case  studies 

(by  Erkan  Giinaydinoglu  (M.Sc.  student  AE  METU,  Asst.  Prof.  Dr.  D.  Funda  Kurtulus  AE  METU) 


In  this  study,  the  flow  structures  and  the  instantaneous  forces  around  SD7003  airfoil 
undergoing  plunging  and  combined  pitching-plunging  motions  are  investigated  numerically 
for  different  frequencies.  The  results  are  a  part  of  the  NATO  AVT-149  Micro  Air  Vehicle 
Unsteady  Aerodynamics  Task  Group  work.  The  resulting  flow-fields  are  visualized  in  terms 
of  normalized  velocity  and  out-of  plane  component  of  vorticity.  The  two-dimensional 
solutions  agree  well  with  experimental  studies.  The  effect  of  Reynolds  number  is  also  studied 
for  a  range  of  10,000  to  60,000.  It  is  observed  that  the  combined  pitching-  plunging  is 
sensitive  to  the  Reynolds  number  while  plunging  is  independent  of  Reynolds  number  in  terms 
of  instantaneous  forces  and  flow  structures. 

INTRODUCTION 


It  has  been  more  than  a  century  that  fluid  mechanists  are  aware  of  that  a  sinusoidal  plunging 
airfoil  is  capable  of  generating  thrust  [18,  19].  Recent  interests  in  Micro  Aerial  Vehicles  take 
into  consideration  the  flapping-wing  motion  as  an  alternative  propulsion  system  for  which  the 
unsteady  low  Reynolds  number  aerodynamics  must  be  studied.  One  of  the  diffuculities  in 
flapping-wing  studies  is  the  transitional  regime  for  Reynolds  numbers  because  of  their  limited 
sizes  as  15  cm  in  height,  length  or  width  [13]  and  the  low-speed  range.  One  other  difficulty  is 
simulating  real  bird-insect  flight  and  that  complex  motion  is  simplified  by  using  basic  models 
i.e.  airfoils-wings  in  pitch,  plunge  or  ramp  motions  [7,  8],  Tuncer  et  al.  [14,  15]  have 
performed  the  Navier-Stokes  computations  to  explore  the  effect  of  flow  separation  on  the 
thrust  and  propulsive  efficiency  of  a  single  flapping  airfoil  in  combined  pitch  and  plunge 
motions.  It  is  observed  that  generated  thrust  and  propulsive  efficiency  are  significantly 
increased  in  case  of  flapping/stationary  airfoil  combination  in  tandem.  Elredge  et  al.  [12] 
analyzed  the  pitch  plunge  motion  with  a  viscous  vortex  particle  method. 


Platzer  et  al.  [10]  emphasizes  that  most  of  the  studies  on  flapping  wing  propulsion  analyzed 
on  literature  are  limited  to  NACA  and  elliptical  airfoils.  The  recent  studies  on  flapping  motion 
are  tent  to  study  on  special  low  Reynolds  number  airfoils  such  as  SD7003  that  is  used  in  the 
study.  Kang  et  a\  [2]  investigated  a  sequence  of  sinusoidal  pitch  and  plunge  motions  of 
SD7003  airfoil  computationally  and  compared  the  results  with  two-dimensional  particle 
image  velocimetry  experiments.  The  forces  are  compared  with  Theodorsen’s  model  and  get 
an  agreement  for  cases  where  massive  separations  do  not  occur.  Yuan  et  al.  [1]  conducted 
Large-eddy  simulations  to  study  plunging  and  pitching-plunging  SD7003  airfoil  and 
compared  the  results  with  PIV  experiments.  01  et  al.  [2]  studied  the  same  motion  numerically 
and  the  results  are  compared  with  PIV  experiments  in  water  tunnel  in  terms  of  vorticity  and 
velocity.  Moreover,  dye  injection  method  is  used  for  a  wide  range  of  Reynolds  number  to 
visualize  the  flow-field  qualitatively.  McGowan  et  al.  [6]  studied  sinusoidal  pure-pitch  and 
pure-p lunge  oscillations  of  SD7003.  The  PIV  experiments  are  compared  with  computations 
with  an  unsteady  RANS  solver  and  Immersed  Boundary  Method  at  Reynolds  number  of 
10,000  and  40,000  which  are  in  the  transitional  region.  In  their  study  experimental- 
computational  agreement  for  plunge  is  satisfied  whereas  the  for  pitch  case  agreement  is  not 
achieved. 

The  flapping  motion  study  has  also  been  performed  experimentally  in  hover  mode  by 
Kurtulus  et  al.  [5]  and  numerically  by  Akay  et  al.  [9].  The  prescribed  flapping  motion  has 
been  investigated  for  different  airfoil  profiles  including  NACA  0012  and  ellipses  with 
different  thicknesses.  Particle  Image  Velocimetry  technique  has  been  used  to  obtain  the  flow- 
field  around  the  NACA0012  airfoil  in  hover  mode  at  Re=1000  and  the  results  are  compared 
with  numerical  solutions  by  Kurtulus  et  al.  [5]. 


In  this  study,  the  unsteady  aerodynamics  of  SD7003  airfoil  undergoing  sinusoidal  plunging 
and  combined  pitching-plunging  motions  is  investigated  at  Reynolds  number  range  of  10,000 
to  60,000  for  reduced  frequencies  of  0.25  and  3.93.  The  numerical  solution  is  compared  with 
the  experimental  data  of  US  Air  Force  Research  Laboratory  (AFRL)  and  TU  Braunschweig 
[2,  8]  and  numerical  solutions  agree  with  experimental  data  in  terms  of  flow  structures  and 
instantaneous  aerodynamic  forces. 


METHOD 


Kinematics 

Pure  plunge  and  combined  pitch-plunge  motions  are  studied  for  two  different  reduced 
frequencies.  The  motions  studied  are  the  case  studies  of  NATO  AVT  149  Micro  Air  Vehicle 
Unsteady  Aerodynamics  Task  Group  .The  equations  of  motion  of  flapping  motion  is  given  in 
equation  (1) 


d(t)  =  Acos(^oot +  ^j  +  a0  (1) 

h(t)  =  h0c  cos(o)t) 

where  h(t )  is  the  plunge  position  of  the  airfoil  with  time  and  6(t)  is  the  pitch  angle  variation 
with  respect  to  the  pivot  point.  The  parameters  for  four  flapping  motion  is  given  in  Table  (1) 
where  k  is  the  reduced  frequency,  k  =  ojc/2Uv,  a0  is  the  mean  angle  of  attack,  h0  is  non- 
dimensional  plunging  amplitude  and  A  is  the  pitch  amplitude  in  degrees. 

Table  1  Parameters  used  to  define  flapping  kinematics 


Cases 

investigated 

k 

a0 

h0 

A 

Re  range 

Pure-plunge 

0.25 

8 

0.5 

0 

10k,  30k,  60k 

Pitch-plunge 

0.25 

8 

0.5 

8.43 

10k,  30k,  60k 

Pure-plunge 

3.93 

4 

0.05 

0 

60k 

Pitch-plunge 

3.93 

4 

0.05 

8.43 

60k 

The  plunging  and  pitching  motion  of  the  airfoil  results  in  an  instantaneous  effective  angle  of 
attack.  Instantaneous  effective  angle  of  attack  for  flapping  motion  is  defined  by  equation  (2) 


aeff 


(  h(t)+  ^ca(t)  cos(a(t)) 
y  Uoo  -  ^ca(t)  sin(a(t)) 


=  a  —  arctan 


(2) 


The  geometric  and  effective  angle  of  attack  variations  of  the  airfoil  during  both  motions  can 
be  observed  in  Figure  1.  The  red  line  denotes  the  effective  angle  of  attack  for  pure-plunge 
motion.  For  pure  plunge  motion,  the  plunge  induced  velocity  and  the  ffee-stream  velocity  are 
used  to  define  the  effective  angle  of  attack.  The  blue  line  in  Figure  1  denotes  effective  angle 
of  attack  with  respect  to  pivot  point  for  pitch-plunge  motion.  The  black  line,  for  which  a  time- 
shift  is  observed  for  both  motions,  is  the  effective  angle  of  attack  for  pitch-plunge  motions 
with  respect  to  the  leading  edge  of  the  airfoil.  In  this  last  case,  the  angular  speed  also  induces 
an  instantaneous  velocity  normal  to  the  chord  direction  at  the  leading  edge  of  the  airfoil. 


Figure  1  The  instantaneous  geometric  and  effective  angle  of  attacks  for  k=0.25,  h=0.5, 
theta=8  cases  (left)  and  k=3.93,  h=0.05,  theta=4  cases  (right) 

Grid  Generation 

In  all  cases  unstructured  O  type  grid  with  outer  diameter  of  60c  around  the  SD7003  airfoil  is 
used.  The  y+  value  for  both  cases  is  chosen  as  less  than  1.  Figure  2  shows  the  whole  domain 
and  the  grid  around  the  airfoil.  A  grid  refinement  process  is  conducted  to  have  grid- 
independent  solution.  The  previous  studies  showed  that  a  grid  with  92358  cells  is  enough  to 
have  the  grid-independent  solution  [3].  Moreover,  960  steps  per  period  are  chosen  for  time 
step. 


Figure  2  Unstructured  grid  distribution  around  SD  7003  airfoil 


Solver  Description 

The  numerical  simulation  of  flapping  motion  is  performed  by  commercially  available  CFD  code 
Fluent  in  which  unsteady,  pressure  based  Reynolds- Averaged  Navier-Stokes  solver  with  second 
order  upwind  spatial  discretization  is  used.  The  pressure-velocity  coupling  is  handled  by  SIMPLE 
algorithm  [17].  The  flapping  motion  of  the  airfoil  and  the  wing  is  implemented  into  the  code  by 
using  ‘dynamic  mesh’  feature  of  Fluent  [11],  To  model  the  motion,  the  whole  grid  is  divided  into 
two  parts:  inner  grid  with  20c  diameter  and  outer  grid  with  60c  diameter.  In  solution  process,  the 
inner  grid  around  the  airfoil  performs  the  prescribed  motion  while  the  outer  grid  is  stationary  and 
deforming  with  appropriate  rules.  The  dynamic  mesh  feature  of  the  code  limits  the  unsteady 
formulation  to  the  first  order  in  time  and  so  for  all  calculations  first-order  temporal  discretization  is 
used.  The  far  boundary  conditions  are  set  as  velocity  inlet  and  pressure  outlet.  The  free-stream  flow 
is  assumed  to  be  turbulent  and  Menter’s  Shear  Stress  Transport  model  is  used  for  turbulence  closure, 
which  is  a  favorable  model  in  prediction  of  adverse  pressure  gradients  [16]. 


RESULTS 

Pitch-plunge  case 

Figure  3  represents  the  u-velocity  component  of  pitch-plunge  with  a  reduced  frequency  of  k=0.25 
and  Re=60,000.  The  results  are  represented  for  different  phases.  The  left  column  represents  the 
numerical  solutions  of  the  present  study.  The  results  are  compared  with  the  TU  Braunschweig  PIV 
solutions  (middle  column)  and  AFRL  solutions  (right  column).  The  PIV  solutions  of  TU 
Braunschweig  is  taken  very  close  to  the  SD  7003  profile,  in  order  to  visualize  the  transition  location 
of  the  airfoil  from  laminar  to  turbulent  flow  with  11  different  PIV  windows  close  to  the  upper 
surface  of  airfoil.  The  numerical  results  seem  to  be  closer  to  their  results  close  to  the  profile. 


Figure  4  shows  the  non-dimensional  vorticity  contours  for  pitch/plunge  case.  At  t/T=0,  the 
instantaneous  angle  of  attack  is  8°.  The  numerical  results  are  very  close  to  the  both  of  the 
experimental  simulations  at  that  instant.  At  t/T=0.25,  the  angle  of  attack  reduces  close  to  zero, 
however  the  effective  angle  of  attack  is  observed  to  be  13.6°.  From  velocity  and  vorticity  contours, 
the  numerical  results  reveal  a  separation  at  this  time  instant  at  the  leading  edge  of  the  airfoil.  Similar 
behavior  is  observed  at  the  experimental  study  in  the  second  column  of  Figure  3  and  Figure  4.  At 
t/T=0.5,  the  angle  of  attack  once  more  increase  to  8°,  for  which  the  effective  angle  of  attack  is  9°.  A 
trailing  edge  vortex  is  observed  at  this  time  instant.  The  leading  edge  vortex  is  separated  from  the 
upper  surface  of  the  airfoil.  At  t/T=0.75,  the  geometric  angle  of  attack  increases  to  16.43°  which  is 
the  maximum  angle  reached  for  this  prescribed  motion.  However,  the  effective  angle  of  attack  at  this 
time  instant  is  close  to  2.4°  so  that  we  observe  an  attached  flow  close  to  the  leading  edge  of  the 
airfoil  at  this  time  instant. 


Figure  3  Normalized  u-velocity  contours  from  computation  ( left  column),  experiment  of  TU 
Braunschweig  (middle  column)  and  experiment  of  AFRL  [2]  (right)  for  pitch/plunge,  Re=60,000, 
k=0.25 
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Figure  4  Normalized  out-of-plane  vorticity  contours  from  computation  (left  column),  experiment  of 
TU  Braunschweig  (middle  column)  and  experiment  of  AFRL  [2]  (right  column)  for  pitch/plunge, 
Re=60,000,  k=0.25 

Pure-plunge  case 

Figure  5  represents  the  u-velocity  component  of  the  pure  plunge  motion  at  k=0.25  and  Re=60,000. 
For  pure  plunge  motion,  the  angle  of  attack  is  8°  however  the  effective  angle  of  attack  changes  as  in 
Figure  1  (left  column).  The  results  are  represented  at  different  phases  of  the  motion  and  they  are 
compared  with  the  experimental  and  numerical  solutions  of  01  et  al  [2],  Pure  plunge  case  represents 
very  good  agreement  with  the  experimental  data.  The  beginning  of  the  motion  at  t=t/T=0  is  very 
close  to  PIV  results  of  01  et  al.[ 2]  where  the  effective  angle  of  attack  is  equal  to  8°.  At  t/T=0.25,  the 
effective  angle  of  attack  increases  to  high  values  on  the  order  of  22°.  This  results  a  massive 
detachment  from  the  airfoil  upper  surface.  U-velocity  component  at  this  time  instant  seems  to  be  in 
very  good  agreement  with  experimental  solutions.  Figure  6  represents  contours  of  normalized  out-of- 
plane  component  of  vorticity  for  pure  plunge  case.  At  t/T=0.25,  a  counter-clockwise  vortex  grows  at 
the  upper  surface  of  the  airfoil,  which  results  the  detachment  of  the  clockwise  leading  edge  vortex 
from  the  airfoil  surface.  At  t/T=0.417,  a  trailing  edge  vortex  in  counter-clockwise  direction  grows, 
and  this  vortex  detaches  from  the  airfoil  surface  very  quickly.  The  trace  of  the  detached  vortex  could 
be  observed  in  Figure  6  at  t/T=0.5.  Similar  trailing  edge  vortex  structure  is  also  observed  from  the 
PIV  measurements  of  01  et  al.  [2], 


10 


Figure  5  Normalized  u-velocity  contours  from  present  computation  (left  column),  experiment  of 
AFRL  [2]  (middle  column)  and  numerical  study  of  Ol  et  al.  [2]  (right  column)  for  pure  plunge, 
Re=60,000  k=0.25 
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Figure  6  Normalized  out-of-plane  component  of  vorticity  contours  from  present  computation  (left 
column),  experiment  of  AFRL  [2]  (middle  column)  and  numerical  study  of  Ol  et  al.  [2]  (right 
column)  for  pure  plunge,  Re=  60,000,  k=0.25 


Reynolds  Number  Effect 

Figure  7  compares  the  u-velocity  components  for  pitch-plunge  case  at  k=0.25  for  different  Reynolds 
numbers  of  10k,  30k  and  60k.  Figure  8  shows  the  contours  of  normalized  out-of-plane  component  of 
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vorticity  for  similar  cases.  At  Re=  10,000,  a  counter-clockwise  vortex  is  observed  at  the  mid- location 
of  the  upper  surface  of  the  airfoil,  which  pushed  the  leading  edge  vortex  towards  up  and  causes  it  to 
detach  from  the  airfoil  surface.  At  t/T=0.333,  when  in  the  same  time  as  this  upper  surface  vortex 
structure,  the  trailing  edge  vortex  in  counter-clockwise  direction  also  grows  for  Re=  10,000  case 
(Figure  8).  The  formation  of  trailing  edge  vortex  is  observed  at  t/T=0.417  for  Re=30,000  and  at 
t/T=0.5  for  Re=60,000. 
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Figure  7  Normalized  u-velocity  contours  for  pitch-plunge  at  Re=l 0,000  (left  column),  Re=3 0,000 
(middle  column)  and  Re=  60,000  (right  column),  k=0.25 
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Figure  8  Contours  of  normalized  out-of-plane  component  of  vorticity  for  pitch-plunge  at  Re=l  0,000 
(left  column),  Re=30,000  (middle  column)  and  Re=  60,000  (right  column),  k=0.25 

The  pure  plunge  case  with  k=0.25  is  also  solved  for  Re=  10,000;  30,000  and  60,000.  The  u- velocity 
contours  for  all  Reynolds  number  variants  are  given  in  Figure  9  and  Figure  10  represents  the 
vorticity  contours  for  the  same  time  instances  of  Figure  9.  The  leading  edge  vortex  structures  are 
much  more  dominant  for  lower  Reynolds  number  case  (Re=  10,000).  As  Reynolds  number  increases, 
the  position  of  the  counter-rotating  vortex  on  the  upper  surface  of  the  airfoil  shifts  towards  the 
leading  edge  of  the  airfoil  (see  t/T=0.25  in  Figure  9).  At  t/T=0.417,  the  core  of  the  trailing  edge 
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vortex  is  observed  to  be  more  far  away  from  the  airfoil  surface  for  Re=  10,000  case.  As  the  Reynolds 

number  increases,  the  trailing  edge  vortex  becomes  closer  to  the  trailing  edge  of  the  airfoil. 

- 1 - 1  I  I  I  I  I  I  I  l  l 

U/IM  0  0  0.1  0.2  0.3  0  4  0.5  0  6  0.7  0  8  0.9  1.0  1.1  1.2  1.3  1  4  1.5 


Figure  9  Normalized  u-velocity  contours  for  pure-plunge  at  Re=l 0,000  (left  column),  Re=3 0,000 
(middle  column)  and  Re=  60,000  (right  column),  k=0.25 
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Figure  10  Contours  of  normalized  out-of-plane  component  of  vorticity  for  pure-plunge  at 
Re=10,000  (left  column),  Re=30,000  ( middle  column)  and  Re=  60,000  (right  column),  k=0.25 


Figure  1 1  represents  the  aerodynamic  force  coefficients  for  pure  plunge  and  pitch/plunge  cases  at 
different  Reynolds  numbers  for  k=0.25.  The  effect  of  Reynolds  number  is  not  abrupt  for  pure  plunge 
case  (right  column  of  Figure  11).  However,  the  peak  location  of  the  pitch/plunge  case  is  found  to 
shift  as  the  Reynolds  number  increases.  A  similar  behavior  is  observed  after  t/T=0.6  for  Re=30,000 
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and  60,000.  Positive  lift  coefficient  values  are  obtained  during  the  whole  period  of  pitch/plunge 
motion. 


Figure  11  Computed  force  coefficients  for  different  Reynolds  numbers,  pitch-plunge  case  (left 
column)  and  pure  plunge  case  (right  column),  k=0.25 


High  Reduced  Frequency  Cases 

The  results  are  also  obtained  for  a  higher  value  of  reduced  frequency  which  is  equal  to  3.93  for  pure 
plunge  and  pitch/plunge  cases.  Figure  12  shows  the  u-velocity  component  of  pure  plunge  case  at 
k=3.93  and  Re=60000.  The  results  are  also  compared  with  the  PIV  solutions  of  01  et  al.  [8],  The 
solutions  at  4  different  phases  of  the  motion  are  represented  in  Figure  12  and  the  results  are  in  good 
agreement  with  experimental  data.  Similar  vortex  structures  are  observed  in  numerical  solutions  of 
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present  study  and  experimental  solutions  of  01  et  al.  [8]  which  could  be  observed  from  the  vorticity 
contours  in  Figure  13.  Figure  14  shows  the  instantaneous  u-velocity  components  of  pitch/plunge  case 
at  k=3.93  and  Re=60000.  As  the  reduced  frequency  increases,  the  vortex  structures  detach  from  the 
trailing  edge  of  the  airfoil  and  seem  to  be  attached  on  the  upper  surface  of  the  airfoil  during  the 
pitch/plunge  case.  A  Karman-vortex  street  is  obtained  from  the  detached  vortex  structure  from  the 
trailing  edge  of  the  airfoil.  At  t=0,  a  clockwise  trailing  edge  vortex  is  observed  which  detaches  from 
the  airfoil  surface  very  quickly  and  leaving  its  place  to  a  counter-clockwise  trailing  edge  vortex  at 
t/T=0.25  (Figure  15).  This  vortex  grows  until  t/T=0.75  and  detaches  from  the  airfoil  surface.  A  new 
clockwise  trailing  edge  vortex  grows  after  that  instant.  During  the  whole  period,  the  flow  is  found  to 
be  attached  to  the  airfoil  upper  and  lower  surfaces.  No  leading  edge  vortex  is  observed  for  this  case. 


Figure  12  Contours  of  normalized  u-velocity  for  pure-plunge  from  computation  (left  column)  and 
experiment  of  Ol  [8]  (right  column)  k=3.93 
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Figure  13  Contours  of  normalized  out-of-plane  component  of  vorticity  for  pure-plunge  from 
computation  (left  column)  and  experiment  of  Ol  [8]  (right  column)  k=3.93 


Figure  14  Contours  of  normalized  u-velocity  for  pitch-plunge  from  computation,  Re=60,000  and 
k=3.93 
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Figure  15  Contours  of  normalized  out-of-plane  component  of  vorticity  for  pitch-plunge  from 
computation,  Re=60,000  and  k=3.93 

Figure  16  shows  the  instantaneous  lift  and  drag  coefficients  for  pure  plunge  and  pitch/plunge  cases  at 
k=3.93  and  Re=60,000.  The  maximum  of  the  lift  coefficient  is  found  to  be  very  high  compared  to 
low  reduced  frequency  case.  Negative  drag  coefficient  values  are  observed  for  pitch/plunge  case 
during  the  whole  period  that  will  result  a  mean  drag  coefficient  of  -0.2677  as  shown  in  Table  2 

Table  2  Computed  mean  force  coefficients  for  investigated  cases 


Cases  investigated 

Mean  Cl 

Mean  Cd 

Pitch-plunge  Re=60k  k=0.25 

0.856 

0.0224 

Pitch-plunge  Re=30k  k=0.25 

0.802 

0.0376 

Pitch-plunge  Re=10k  k=0.25 

0.683 

0.0609 

Pure-plunge  Re=60k  k=0.25 

0.766 

0.0805 

Pure-plunge  Re=30k  k=0.25 

0.779 

0.0931 

Pure-plunge  Re=10k  k=0.25 

0.682 

0.0956 

Pitch-plunge  Re=60k  k=3.93 

0.838 

-0.2677 

Pure-plunge  Re=60k  k=3.93 

0.703 

-0.0664 
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Figure  16  Computed  force  coefficients  for  high  frequency  case,  Re=60,000  and  k=3.93 


The  pitch  plunge  and  pure  plunge  motions  are  also  investigated  with  a  flat  plate  having  2.5% 
thickness.  The  normalized  u-velocity  component  and  vorticity  contours  at  different  phases  of  the 
pitch/plunge  motion  are  given  in  Figure  17  at  Re=60000.  Figure  18  shows  the  normalized  u-velocity 
contours  and  vorticity  contours  for  pure  plunge  case  at  different  phases  of  the  motion. 
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Figure  17  Contours  of  normalized  u-velocity(left)  and  vorticity  (right)  for  pitch/plunge  motion  of  a  flat  plate  with 
2.5%  thickness 
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tp=270° 

Figure  18  Contours  of  normalized  u-velocity(left)  and  vorticity  (right)  for  pure  plunge  motion  of  a  flat  plate  with 


tp=150° 


2.5%  thickness 
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CONCLUSIONS 


The  pitch  plunge  and  pure  plunge  cases  at  different  reduced  frequencies  are  investigated  during  this 

study.  The  results  are  compared  with  the  experimental  solutions  available  in  literature.  The  solutions 

are  part  of  the  AVT-149  technical  team  work. 
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1.2  Steady  state  computations  of  SD  7003  and  flat  plate  with 
3.5%  thickness 


a)  Steady  state  solution  around  SD  7003  at  Re=  60000 


Figure  17  and  Figure  18  show  the  grid  used  for  the  stead  state  solution  around  SD  7003  airfoil. 


Figure  19  Unstructured  grid  around  SD  7003 


Figure  20  Boundary  layer  close  to  the  airfoil 
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Figure  21  and  Figure  22  show  the  lift  and  drag  coefficients  for  SD  7003  airfoil  at  Re=60000.  The 
numerical  results  of  the  present  study  are  compared  with  the  experimental  and  numerical  solutions 
in  the  literature. 


Figure  21  Lift  coefficient  versus  angle  of  attack  at  Re=60000  for  SD  7003  airfoil  [compared  with  from  Shyy,  W., 
Lian,  Y.,  Tang,  J.,  Viieru,  D.  and  Liu,  H.  (2008)  Aerodynamics  of  Low  Reynolds  Number.  Flyers,  Cambridge 
University  Press,  New  York] 


Figure  22  Drag  coefficient  versus  angle  of  attack  at  Re=60000  for  SD  7003  airfoil 
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Figure  23  shows  the  pressure  coefficient  contours  around  SD  7003  airfoil  at  Re=60000  for  different 
angles  of  attack.  As  the  angle  of  attack  increses,  the  suction  peak  region  close  to  the  leading  edge 
of  the  airfoil  spreads  towards  the  trailing  edge  of  the  airfoil.  Stall  angle  is  close  to  1 1  degrees.  Figure 
24  shows  the  non-dimensional  u-velocity  component  around  SD  7003  at  different  angles  of  attack. 
The  non-dimensional  vorticity  contours  are  also  represented  in  Figure  25  for  the  stead-state 
solution. 


29 


Figure  23  Pressure  coefficient  ( Cp)  around  SD  7003  airfoil  at  different  angle  of  attacks. 
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Figure  24  (U  /  U^)  graph  around  SD  7003  at  different  angles  of  attack  (Re=60000). 
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a  =  16 
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Figure  25  Non  dimensional  vorticity  (  COzC  /  U^)  graph  around  SD  7003  airfoil. 
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Figure  26  Pressure  coefficient  around  SD  7003  airfoil  at  different  angles  of  attack 


Figure  27  represents  the  streamlines  around  SD  7003  airfoil.  A  trailing  edge  vortex  structure  is 
obseved  at  12°  angle  of  attack.  A  leading  edge  vortex  structure  detached  form  the  airfoil  surface 
with  a  trailing  edge  vortex  is  observed  at  16°  angle  of  attack.  The  suction  region  on  the  upper 
surface  of  the  airfoil  becomes  dominant  after  10°  angle  of  attack  as  can  be  observed  from  Figure 
24. 


Figure  27  Streamlines  around  SD  7003  airfoil  for  12°  and  16°  angles  of  attack  (Re=60000). 
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b)  Steady  state  solution  around  Hat  plate  with  3.5%  thickness  at  Re=  60000 


For  steady  state  solutions  around  flat  plate,  84142  elements  are  used  for  the  grid.  The  lift  and 
drag  coefficients  at  different  angles  of  attack  are  given  in  Table  3  for  Re=60000. 

Table  3  Steady-state  lift  and  drag  coefficients  for  different  angles  of  attack  at  Re=60000 


AOA 

CL 

CD 

4 

0.46432 

0.044677 

8 

0.78616 

0.12399 

12 

0.68309 

0.16936 

16 

0.85031 

0.26908 

u/U^:  0  0.1  0.2  0.3  0.4  0.5  0.6  0  7  0.8  0.9  1  1.1  1.2  1.3  1.4  1.5 
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or  =  16 
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0  0.5  1  1.5  2  ’n  0  0.5  1  1.5 
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Figure  28  u/U^  contours  around  flat  plate  with  3.5%  thickness  (Re=60000) 


34 


007cAJ  ■  -36  -30  -24  -18  -12  -6  0  6  12  18  24  30  36 

Z  CO 


Figure  29  Non-dimensional  vorticity  contours  around  flat  plate  with  3.5%  thickness  (Re=60000) 
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Figure  30  Pressure  coefficient  contours  around  flat  plate  with  3.5%  thickness  (Re=60000) 


Figure  31  Pressure  coefficient  contours  on  the  surface  of  flat  plate  with  3.5%  thickness  (Re=60000) 
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1.3  Unsteady  Aerodynamic  Computations  of  Different  Sinusoidal 
Flapping  Motions 

(by  Busra  Akay  (Ph.D.  student  TU  Delft,  Asst.  Prof.  Dr.  D.  Funda  Kurtulus  AE  METU) 

This  study  aims  to  analyze  the  present  numerical  method  by  simulating  two  dimensional  flapping 
motion  in  hovering  mode.  Direct  Numerical  Simulation  is  used  to  solve  the  flow  field  around  the  two 
dimensional  wing  during  the  flapping  motion.  Unsteady,  laminar,  incompressible  two  dimensional 
Navier-Stokes  Equations  are  solved  by  using  moving  grid  technique.  Two  different  sinusoidal 
motions  are  simulated  by  implementing  the  sinusoidal  translational  and  angular  motions  using  the 
same  numerical  solver  and  compared  with  the  experimental  results.  This  study  is  performed  to 
explore  the  computational  performance  of  the  present  numerical  method  and  to  analyze  the 
sinusoidal  flapping  motion  aerodynamics  for  different  Reynolds  numbers  in  the  range  lOMO3  by 
using  different  kinematics. 

INTRODUCTION 

Interest  in  the  aerodynamics  of  insect  flights  has  increased  in  conjunction  with  the  concept  of  Micro 
Air  Vehicles  (pAVs).  Based  on  their  size,  flying  insects  operate  in  a  wide  range  of  Reynolds 
numbers;  from  approximately  10  to  105  [1].  Operating  Reynolds  numbers  range  of  pAVs  is  similar 
with  those  of  birds  or  insects.  Thus,  this  similarity  leaded  most  of  the  researchers  to  understand  the 
aerodynamic  basis  for  the  flight  of  birds  and  insects. 

Wang  has  defined  hovering  as  an  extreme  mode  of  flight  where  the  forward  velocity  is  zero.  To 
hover  insects  must  draw  clean  air  from  the  ambient  flow  and  get  rid  of  the  'messy  vortices'  they 
have  created  to  obtain  a  large  periodic  lift  [2].  However,  the  hovering  flight  is  quite  expensive.  While 
weak  fliers  and  strong  flying  birds  invest  about  15  and  20  percent  (respectively)  of  the  total  body 
weight  in  the  breast  muscles,  hummingbirds  invest  about  30  percent  of  the  total  body  weight  in  the 
breast  muscles  [3],  Although  flapping  wing  design  is  more  complex  than  a  fixed  wing  design,  there 
are  many  reasons  to  explore  the  possibilities  of  flapping  wing  flight.  While  the  vehicle  becomes 
smaller,  the  fixed  wing  application  becomes  less  reasonable.  The  lift  which  a  fixed  wing  generates 
to  support  the  weight  of  the  vehicle  is  directly  proportional  to  wing  area  and  velocity  of  air  flow  over 
the  wing.  Thus,  the  smaller  the  vehicle,  the  less  lift  it  can  supply.  Provided  with  enough  power,  a 
vehicle  with  flapping  wings  could  actually  takeoff  and  land  vertically  [4]. 

To  explore  this  complex  topic,  many  researches  have  been  performed  numerically  and 
experimentally.  Here,  some  of  the  numerical  works  are  reviewed.  Many  applications  have  been 
performed  by  using  Navier-Stokes  solvers  in  the  field  of  flapping  motion.  NACA  0012  airfoil  profile 
has  been  used  as  wing  section  in  many  applications  ([5]-[9]).  Elliptic  profiles  have  also  been  used  to 
investigate  the  flapping  motion  characteristics  ([10]-[13]).  Lan  and  Sun  [10]  explored  the  flapping 
motion  at  Re=1000  by  using  a  Navier-Stokes  solver  for  incompressible  flow  implementing  moving 
overset  grid.  The  results  show  that,  if  the  insect  employs  a  larger  angle  of  attack  or  changes  the 
timing  of  wing  rotation,  much  greater  lift  can  be  produced  for  maneuvering  and  for  other  purposes. 
3D  flapping  motion  of  the  model  fruit  fly  wing  at  Re=136  has  been  investigated  by  Sun  and  Tang 
[11]  with  some  insights  into  the  unsteady  aerodynamic  force  generation  process  from  the  force  and 
flow-structure  information.  They  compared  their  results  with  the  model  wing  experimental  results 
and  fruit-fly  data  provided  by  Dickinson  et  al.  [14]  and  Weis-Fogh  [15].  Weis-Fogh  [15]  aimed  to 
provide  new  material  and  novel  solutions  to  make  use  of  the  large  number  of  observations  on  freely 
flying  animals.  His  major  conclusion  is  that  most  insects  perform  normal  hovering  on  the  basis  of  the 
well-established  principles  of  steady-state  flow.  However,  one  must  also  realize  that  any  type  of 
flapping  flight  involves  also  non-steady  periods,  particularly  at  the  reversal  points  where  active 
pronation  and  supination  occur.  Wang  ([2],  [12])  has  analyzed  2D  hovering  and  flapping  flight  on 
elliptic  wing  section  to  identify  the  vortex  shedding  and  their  frequencies  in  the  Re  number  range  of 
102<Re<104.  Eldredge  [13]  has  performed  DNS  solutions  with  viscous  vortex  particle  method  to 
investigate  the  pitching  and  plunging  motion  at  Re=550. 
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NUMERICAL  METHOD 

Direct  Numerical  Simulation  (DNS)  technique  is  used  to  solve  the  present  flapping  motion  problems. 
Simulation  is  performed  for  laminar,  incompressible  flow  condition.  The  CFD  code  (Star-CD)  used  in 
this  simulation  has  the  capability  to  solve  transient  flow  problems,  to  use  moving  mesh  with  arbitrary 
motions,  to  handle  user  defined  properties  and  conditions  by  the  use  of  user-defined  subroutines, 
and  it  also  has  the  capability  of  handling  a  large  variety  of  boundary  conditions,  and  offers  a  range 
of  moving  mesh  features. 

Wing  Models  and  Their  Kinematics 

Two  different  sinusoidal  flapping  motions  are  analyzed  in  the  present  study.  The  first  one  is  defined 
by  Wang  et  al.  [16]  and  in  this  study  it  is  called  Type  A’  flapping  motion.  It  is  applied  to  an  ellipse 
having  12%  thickness  and  0.01m  chord  length.  The  computed  results  are  then  compared  with  three 
dimensional  experiments  and  empirical  data  [16].  The  second  flapping  motion  is  defined  by 
Freymuth  [17]  and  in  this  study  it  is  called  Type  B’  flapping  motion.  It  is  implemented  to  an  elliptic 
wing  model  having  a  thickness  of  1.6mm  and  a  chord  length  of  0.0254m  as  in  Freymuth’s 
experiments  [17].  Defined  boundary  conditions  and  inner  grid  domains  used  in  the  both  simulations 
are  presented  in  Figure  32.  Grid  domain  is  formed  via  GRIDGEN  VI 5,  a  package  programmed  to 
generate  grid  domain.  O-type  grid  domain  is  used  around  the  profiles.  According  to  the  results  of 
grid  refinement  study,  229x340  (229  number  of  nodes  around  the  profile)  grid  domain  is  used  in  the 
numerical  simulation  of  this  study.  The  domains  are  formed  in  two  sub-domains,  inner  domain  is 
finer.  The  radius  of  the  whole  domain  is  20  chords  having  a  total  of  77292  cells  for  the  case  of  Type 
A’  simulation.  259  numbers  of  nodes  are  put  on  the  profile  of  case  of  Type  B  simulation.  The  radius 
of  the  whole  domain  is  20  chords  having  82560  cells. 


Symmetry 

Boundaries 


Figure  32  Defined  boundary  conditions  and  inner  grid  domains.  Profile  of  (a)  Type  A  simulation  (b)  Type  B 
simulation 
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Kinematics  of  ‘Type  A’  Flapping  Motion 

For  the  simulation  of  flapping  kinematics  of  Wang  et  al.  [16],  an  ellipse  of  12%  chord  thickness  is 
used  (c=0.01m).  For  the  ellipse  (e=12%c),  the  same  grid  domain  is  used  as  the  previous  study.  The 
wing  follows  a  sinusoidal  flapping  and  pitching  motion  (Eqs.  1-2,  respectively)  [16].  Specifically,  the 
wing  sweeps  in  the  horizontal  plane  and  pitches  about  its  spanwise  axis  with  a  single  frequency  f: 


x(t)  = 


"~~cos(2^z/r) 


(1) 


a(t)  =  a0  +  P  sin(2^  +  tf>) 


(2) 


where  x(t)  is  the  position  of  the  center  of  the  wing,  and  a  (t)  is  the  wing  orientation  with  respect  to 
the  x-axis.  By  definition,  the  translational  and  angular  velocities  are  given  by  U0(t)=dx(t)/dt  and  Q 
(t)=da(t)/dt.  The  parameters  include  the  stroke  amplitude  A0,  the  initial  angle  of  attack  a0,  the 
amplitude  of  pitching  angle  of  attack  /?,  the  frequency  f  and  the  phase  difference  <j>  between  x(t)  and 

a(t) . 


Upstroke 


Dmvnstroke 


Figure  33  Sinusoidal  motion  of  the  profile  during  one  stroke. 

The  translational  motion  of  the  wing  is  completely  specified  by  two  dimensionless  parameters, 
Reynolds  number,  Re=\Jmaxcl  v=  7tfA0c /  v ,  and  A0lc,  where  Umax  is  the  maximum  flapping  velocity, 

and  c  is  the  chord  length.  From  their  steady-state  2D  numerical  data  Wang  et  al.  [16]  found  the 
approximated  empirical  correlations  for  both  of  the  aerodynamic  coefficients,  namely,  C/.(Eq.  3)  and 
CD  (Eq.  4)  in  terms  of  angle  of  attack  a . 


C,  =  1.2sin(2  a) 

(0) 

CD  =1.4-  cos(2cr) 

(4) 

The  constants  depend  on  the  Reynolds  number,  details  of  the  wing,  etc.  They  implemented  this 
empirical  data  (Eq.  3-Eq.  4)  for  all  of  the  instantaneous  angle  of  attack  variations  that  they  have 
investigated.  For  each  a  value,  CL  and  CD  values  are  calculated.  Quasi-steady  translational  lift  ( LT ) 
and  drag  forces  ( DT )  are  calculated  0.5pu2CL  and  0.5pu2CD,  respectively.  All  of  the  numerical  and 
empirical  forces  are  normalized  by  the  maxima  of  the  corresponding  to  quasi-steady  forces  as 
described  in  the  study  of  Wang  et  al.  [16]. 
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Kinematics  of  ‘Type  B’  Flapping  Motion 

In  this  part  of  the  study,  the  flapping  motion  as  defined  by  Freymuth  [17]  is  investigated.  Freymuth 
[17]  used  a  planar  airfoil  having  a  thickness  of  1.6mm  and  a  chord  of  c=2.54cm  with  rounded  edges 
to  execute  the  combined  plunging  and  pitching  motions  in  the  experiments.  In  the  present  numerical 
investigations  an  elliptical  wing  having  the  same  thickness  and  chord  as  Freymuth’s  model  is  used. 

The  airfoil  performs  a  translating  (plunging)  motion  [17]  h  in  horizontal  direction  (Eq.  5): 


h  =  ha  sin(2  rfi) 


(5) 


where  ha  is  the  amplitude  of  linear  translation,  f  is  the  frequency  of  sinusoidal  oscillation  and  t  is  the 
time.  Considering  that  the  airfoil  performs  a  pitching  motion  (Eq.  6)  simultaneously,  around  its  half 
chord  axis: 


a  =  a0  +  aa  sm{2nft  +  tf>) 


(6) 


where  a  is  the  pitch  angle  with  respect  to  the  horizontal  as  shown  in  Figure  34,  a0  is  the  mean  pitch 
angle,  aa  is  the  pitch  amplitude  and  <j>  is  the  phase  difference  between  pitching  and  plunging. 

Dimensionless  parameters  of  the  system  are:  a0,  the  dimensionless  plunge  amplitude  hjc  and 
a  Reynolds  number; 


Rf  =  2  7tfhadv 


(7) 


based  on  the  maximum  plunge  speed  2;z/7tfland  on  c,  where  v  is  the  kinematic  viscosity. 

Two  simple  modes  of  hovering  were  initially  identified  by  Freymuth  [17],  In  the  present  study  only 
Mode  1  has  been  analyzed.  “Mode  1”  or  “water  treading  mode”  is  characterized  by  a0=0°  and  </>=  90° 

and  is  sketched  in  Figure  34.  The  airfoil  starts  a  cycle  from  the  position  of  having  pitch  amplitude 
(aa)  at  middle  of  the  downstroke  (indicated  as  right  arrow).  It  moves  a  distance  2ha  to  the  right  to 
reach  its  initial  position.  The  right  edge  of  the  airfoil  is  leading  during  its  motion  to  the  right  but  when 
the  airfoil  returns  left,  the  leading  and  trailing  edges  switch  their  roles. 


-e 


Figure  34  Sketch  of  combined  translating-pitching  motions  of  the  airfoil  for  one  cycle  of  mode  1  hovering  (ao=0°, 
^=90°). 


To  characterize  the  time  averaged  thrust  T  on  the  airfoil,  a  thrust  coefficient  CT  is  defined  in 
Reference  [1 7]  (Eq.  8): 
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T 


(8) 


CT 


0.5  pVt2cl 


where  p  is  the  air  density,  l»c,  is  the  span  of  the  airfoil. 

RESULTS 

This  part  is  devoted  to  the  numerical  analysis  of  different  flapping  motion  kinematics  for  different 
Reynolds  numbers  in  the  range  of  (lOMo3).  To  explore  the  consistency  of  the  present  numerical 
method,  experimental  results  obtained  from  the  previous  studies  in  the  literature  are  compared  with 
the  numerical  results  of  the  present  study. 

Analysis  of  Type  A  Flapping  Motion 

Sinusoidal  angle  of  attack  and  velocity  distributions  of  the  motion  for  different  stroke  amplitudes  are 
presented  in  Figure  35  to  show  the  difference  of  the  kinematics  by  the  change  of  A0/c  value.  The 
frequency  of  the  former  is  approximately  twice  of  the  later  one. 


Figure  35  Instantaneous  velocity  and  angle  of  attack  (  CC  (t))  distribution  vs.  time 


In  the  following  results,  time  is  non-dimensionalized  with  respect  to  the  flapping  period  of  the  case. 
In  Figure  36,  the  computed  forces  obtained  from  the  present  study  are  compared  with  the 
experimental  and  empirical  data  of  Wang  et  al.  [16].  The  empirical  data  is  carried  out  by  using  the 
Eqs.  3-4  based  on  the  translational  velocity.  The  forces  are  normalized  by  the  maxima  of  the 
corresponding  quasi-steady  forces  [16],  For  lift  coefficient  distribution,  although  there  is  a  slight  over 
estimation  at  the  mid-strokes,  the  present  computations  catch  the  CL  values  of  the  experiment 
during  translation  of  the  wing  (see  Figure  36).  Also,  the  present  computations  estimate  the  right 
rotation  position  as  in  the  experiment.  For  drag  coefficient,  it  is  observed  that  the  general  behavior 
of  the  distribution  obtained  by  the  present  computations  is  very  similar  with  that  of  the  experiment 
[16].  Especially  in  the  translational  phase,  they  are  very  successful.  Generally,  it  is  noted  that  the 
present  computations  are  very  good  at  estimating  the  force  coefficients  of  this  problem. 
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Figure  36  Lift  and  Drag  coefficients  comparison  between  experiment  [16],  present  computation  and  quasi-steady 
estimations  [16]  for  symmetric  ((p=0)  rotation,  Re=75,  and  A0/c=2.8.  Time  is  non-dimensionalized  with  the 
flapping  period  of  the  case  1. 

A  comparison  between  computation  and  experiments  of  Wang  et  al.  [16]  and  the  present 
computational  results  will  be  presented  for  the  case  of  A0/c=4.8  and  Re=115.  In  the  first  column  of 
Figure  37,  2D  computational  results  [16];  in  the  second  column  3D  experimental  results  in  a  2D  slice 
at  0.65R  taken  from  DPIV  measurements  [16]  and  in  the  third  column  2D  computational  results  of 
the  present  study  are  represented  respectively  in  columns  1  to  3  of  Figure  37.  Ten  different  time 
sequences  are  shown  during  the  fourth  stroke  for  each  case.  Time  is  non-dimensionalized  with  the 
flapping  period  of  the  case. 


The  vorticity  contours  are  presented  (see  Figure  37)  to  show  the  major  features  of  vortex  dynamics 
through  a  complete  stroke  cycle.  When  the  three  results  at  each  indicated  time  steps  are  compared, 
it  is  seen  that  the  major  features  of  vortex  dynamics  are  similar.  The  color  scale  for  vorticity  of 
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computation  and  experiments  [16]  did  not  correspond  to  the  exact  same  contour  values.  Due  to  lack 
of  the  color  scale  of  vorticity  contours  achieved  by  Wang  et  al.  [16],  the  figure  of  our  computational 
results  and  their  results  should  be  viewed  more  qualitatively  than  quantitatively. 

Notice  that  even  though  the  kinematics  of  left  and  right  strokes  is  identical,  the  flow  fields  for  each 
case  differ  slightly  for  the  computational  results  of  Wang  et  al.  [16],  However,  this  discrepancy  in  the 
flow  field  is  not  observed  in  our  computational  results. 


Figure  37  Instantaneous  vorticity  contour  for  case  of  A0/c=  4.8,  Re=115,  <p=0.  First  two  columns  are  the  results  of 
Wang  et  al.  [16]  and  the  third  column  is  the  results  of  the  present  study. 
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6  (Cont’d)  Instantaneous  vorticity  contour  for  case  of  A0/c=  4.8,  Re=115,  (p=0.  First  two  columns  are  the 


results  of  Wang  et  al.  [16]  and  the  third  column  is  the  results  of  the  present  study. 


Analysis  of  Type  B  Flapping  Motion 

Sinusoidal  variations  angle  of  attack  and  velocity  distributions  of  the  motion  for  mode  1  are 
presented  in  Figure  38.  Mode  1  flapping  motion  prescribed  by  Freymuth  [17]  is  applied  to  the  wing 
model.  In  Figure  39,  force  coefficient  computed  by  using  the  present  numerical  method  is  compared 
with  the  experimental  data  of  Freymuth  [17].  The  computed  lift  is  non-dimensionalized  according  to 
Eq.  8.  The  results  should  be  compared  by  bearing  in  mind  the  differences  between  two  studies  such 
as  wing  model,  three  dimensional  effects  of  the  experiments  [17], 
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Figure  38  Instantaneous  angle  of  attack  (a(t))  and  velocity  distribution  vs.  time  for  mode  1  hovering. 
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Figure  39  Lift  coefficient  comparison  between  experiment  [17]  and  present  computation  for  mode  1,  aa=66°,  a0=0, 
(p=90  and  R^1700. 

Entire  cycle  of  airfoil  flapping  in  Mode  1  with  aa=66°,  ha/c=1.5  is  shown  in  Figure  40  (a-b).  Black  and 
white  pictures  belong  to  Freymuth’s  experiment  [17].  Flow  visualization  was  by  means  of  the 
titanium-tetra-chloride  method  described  by  Freymuth  et  al.  [18].  Colorful  pictures  are  results  of  the 
present  study.  Frames  are  ordered  into  columns  from  top  to  bottom  and  columns  are  ordered  from 
left  to  right.  Time  between  consecutive  frames  is  At=1/16s.  One  should  analyze  the  pictures  by 
following  the  first  and  second  columns  together.  First  row  of  left  two  columns  show  the  farthest  right 
position  of  the  airfoil  (lower  right  corner  of  the  frames).  From  this  position  to  the  bottom  of  the  third 
and  fourth  columns,  the  airfoil  moves  its  left  position.  During  this  movement,  airfoil  creates  a 
clockwise  (blue)  rotating  vortex  which  is  very  similar  with  the  experiment  [17]  (see  Figure  40  (a)- 
indicated  in  the  last  row).  In  columns  1  and  2  of  Figure  40  (b),  the  airfoil  moves  to  the  right.  The 
previously  generated  clockwise  (blue)  rotating  vortex  (see  Figure  40  (a)-indicated  in  the  last  row) 
starts  to  detach  from  the  upper  surface  of  the  airfoil  and  moves  upward.  And  a  new 
counterclockwise  (red)  rotating  vortex  is  formed  and  grows  (see  Figure  40  (b)-indicated  in  the  last 
row).  This  process  repeats  during  each  cycle  and  results  in  an  upward  moving  vortex  street.  These 
vortex  formations  were  observed  in  the  experiment  also  [17]. 
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Figure  40  (a)  Instantaneous  vorticity  contours  for  Mode  1.  Black  and  white  pictures  are  results  of  Freymuth  [17] 


and  the  colorful  pictures  are  results  of  the  present  study.  aa=66°,  ha/c=1.5,  Rf=340,  f=lHz.,  At=l/16s. 
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Figure  40  (b)  Instantaneous  vorticity  contours  for  Mode  1.  Black  and  white  pictures  are  results  of  Freymuth  [17] 
and  the  colorful  pictures  are  results  of  the  present  study.  aa=66°,  ha/c=1.5,  Rf=340,  f=lHz.,  At=l/16s. 


CONCLUSION 

A  study  is  performed  to  analyze  the  different  sinusoidal  flapping  motion  kinematics  defined  by  Wang 
et  al.  [16]  and  Freymuth  [17]  for  different  Reynolds  numbers  in  the  range  of  lOMo3  by 
implementation  of  the  sinusoidal  translational  and  angular  motions.  Force  coefficients  and  vorticity 
contours  obtained  from  the  experiments  [16],  [17]  and  present  study  are  compared. 


47 


The  sinusoidal  motion  defined  by  Wang  et  al.  [16]  is  applied  to  a  thin  wing  element  of  elliptic  cross 
section  having  a  thickness  of  12%  of  its  chord.  The  computed  forces  and  vorticity  distribution 
obtained  from  the  present  study  are  compared  with  3D  experimental,  2D  numerical  and  empirical 
data  of  Wang  et  al.  [16].  It  is  observed  that  the  present  computational  method  is  good  at  estimating 
the  force  coefficients  and  the  major  features  of  vortex  dynamics  of  the  problem. 

The  sinusoidal  flapping  motion  defined  by  Freymuth  [17]  is  implemented  to  an  elliptic  profile  having 
1.6mm  thickness  and  0.0254m  chord  length.  Force  coefficients  and  vortex  dynamics  obtained  from 
the  experiments  of  Freymuth  [17]  and  present  study  are  compared.  Although  some  parameters  are 
different  for  numerical  and  experimental  tests  good  agreement  is  observed  between  these  two 
studies. 
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2.  Experimental  Method 


2.1  Vortex  topology  using  PIV 

(by  Dr.  D.  Funda  Kurtulus  AE  METU) 


The  flow  around  a  flapping  airfoil  is  investigated  for  different  angles  of  attack  by  means  of  particle  image 
velocimetry  (PIV)  technique  at  Re=1000.  The  objective  is  to  show  the  influence  of  flapping  motion  in  a 
flow  where  there  exist  abundant  vortex  structures.  The  instantaneous  circulations  are  calculated  for  the 
vortices  generated  at  the  leading  edge  and  trailing  edge  of  the  airfoil  during  the  impulsive  start  of  the 
flapping  motion  and  during  the  pseudo-periodic  motion  where  the  influence  of  the  impulsive  start 
disappears.  The  vortex  generation  mechanism  at  the  impulsive  start  is  studied  for  the  understanding  of  the 
flapping  motion  vortex  topology  in  order  to  distinguish  the  characteristics  of  vortices  generated  during  the 
motion  without  the  influence  of  the  traces  from  the  previous  strokes.  By  comparison  with  the  same 
location  for  the  periodic  motion,  the  influences  of  the  traces  of  the  previous  motion  and  the  last  vortices 
shedded  at  the  end  of  the  translational  phase  are  exhibited. 


1.  Introduction 

In  order  to  understand  the  high  lift  mechanisms  of  flapping  wings,  lots  of  investigations  are  performed  these 
last  years.  The  vortical  flow  patterns  in  the  wake  of  a  NACA  0012  airfoil  pitching  at  small  amplitudes  are 
studied  experimentally  by  Koochesfahani  (1989)  in  a  low-speed  water  channel  by  considering  the  effect  of 
both  sinusoidal  and  non-sinusoidal  shape  of  the  waveform.  The  unsteady  flow  structures  are  investigated  with 
dye  visualization  by  Freymuth  (1990).  He  has  worked  on  experimental  simulation  of  sinusoidally  pitching  and 
plunging  thin  flat  plate  to  investigate  the  thrust  generation  mechanism  in  hover  mode  with  Re  number  range  of 
340  to  1700.  He  examines  the  dynamic  stall  vortices  effects  in  order  to  obtain  large  thrust  values. 

PIV  and  SPIV  have  been  often  used  as  an  experimental  tool  to  analyze  the  flapping  flight  vortical 
structures.  Poelma  et  al.  (2006)  have  performed  a  3D  Steoreoscopic  PIV  experiment  in  a  mineral  oil 
tank  for  the  flapping  flight  with  dynamically  scaled  robotic  wing  at  Re=256.  Wang  et  al.  (2004)  have 
performed  PIV  experiment  together  with  the  force  measurement  by  using  2D  force  sensors  at  low  Re 
number  (50<Re<200).  Dynamically  scaled  robotic  fly  has  been  used  in  their  experiments.  Laser 
sheet  visualizations  and  2D  PIV  experiments  in  a  water  tank  have  been  performed  by  Kurtulus  et  al. 
(2008)  with  NACA  0012  airfoil  in  the  Re  number  range  of  500<Re<2000.  01  (2007)  and  Radespiel 
et  al.  (2007)  recently  investigated  vortical  structures  generated  during  the  pitch  and  plunge  motions 
of  SD7003  airfoil  at  Re=60000. 

Dickinson  et  al.  (1999)  have  performed  an  experiment  in  a  mineral  oil  tank  with  a  dynamically 
scaled  model  of  the  fruit  fly  to  investigate  3D  flapping  motion  at  Re=136.  The  wing  was  equipped 
with  a  2D  force  transducer.  Before  this  study,  Dickinson  and  Gotz  (1993)  have  performed  similar 
experiment  in  an  aquarium  by  using  2D  impulsively  moved  model  wings  in  the  range  of 
10<Re<1000.  In  this  experiment,  the  direct  force  measurements  were  correlated  with  CCD  camera 
visualization.  Clap  and  fling  movement  have  been  analyzed  using  dynamically  scaled  mechanical 
model  of  the  small  fruit  fly  by  Lehmann  et  al.  (2005).  They  performed  3C  SPIV  experiment  and  used 
force  transducer  to  get  the  force  distribution  in  the  range  of  100<Re<200. 

Some  zoological  experiments  have  also  been  carried  out  by  the  researchers.  Tian  et  al.  (2006)  have 
implemented  PIV  by  using  fog  generator  in  a  flight  cage  to  understand  3D  kinematic  motion  during 
straight  and  turning  of  bat  flight  in  the  Re  number  range  of  104<Re<105.  Galvao  et  al.  (2006)  have 
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explored  3D  mammalian  flight  with  compliant  membrane  wing  models  in  the  range  of 
70000<Re<200000.  They  have  used  low-speed,  low  turbulence  wind  tunnel  with  stereo 
photogrammetric  method.  2D  biomimetic  flapping-pitching  wing  is  analyzed  by  Singh  et  al.  (2005) 
using  laser  sheet  visualization  method.  Usherwood  et  al.  (2005)  have  investigated  flapping  flight  to 
obtain  dynamic  pressure  maps  of  wings  and  tails  of  Pigeons.  Platzer  and  Jones  (2006)  have  recently 
reviewed  different  flapping  wing  studies. 

In  the  view  of  all  these  different  studies,  different  effects  are  analyzed  in  order  to  explain  the  fluid  physics 
responsible  for  lift  generation  (Dickinson  1994,  Sane  2003,  Viiem  et  al.  2006).  Wagner  effect  (Wagner  1925, 
Sane  2003)  is  only  applicable  at  the  beginning  of  the  impulsive  start  of  the  motion.  When  an  inclined  wing 
starts  impulsively  from  rest,  the  circulation  around  it  does  not  immediately  attain  its  steady  state  value.  This 
effect  is  found  to  be  not  strong  both  in  2D  experiments  of  Dickinson  and  Gotz  (1993)  and  3D  flapping 
translation  experiments  of  Dickinson  et  al.  (1999)  at  low  Re  number  regimes. 

Delayed  stall  effect  results  from  the  growth  of  the  attached  leading  edge  vortex  which  produces  very  high  lift 
coefficients  for  several  chord  length  translations  before  its  detachment  from  the  airfoil  upper  surface  and 
before  the  formation  of  trailing  edge  vortex.  These  leading-edge  vortices  result  in  a  large  suction  region  on  the 
upper  surface  of  the  airfoil  during  translation  phase  of  the  flapping  motion  and  enhance  the  lift  (Kurtulus  et  al. 
2008).  This  phenomenon  repeats  by  the  detachment  of  trailing  edge  vortex  and  creation  of  new  leading  edge 
vortex  (von  Karman  vortex  street). 

Kramer  effect  (rotational  forces  or  fast  pitching  up)  is  due  to  the  rotation  of  the  airfoil  at  the  end  of  the  stroke 
reversals.  The  wing  generates  a  rotational  circulation  (proportional  to  the  angular  velocity  of  rotation)  in  the 
fluid  to  counteract  the  effects  of  rotation  and  as  a  result  rotational  force  peaks  are  generated  during  stroke 
reversals  to  re-establish  the  Kutta  condition  during  the  rotation.  (Sun  and  Tang  2002,  Sane  and  Dickinson, 
2002,  Sane  2003,  Viiem  et  al.  2006).  Another  phenomena  is  the  added  mass  effect  due  to  the  acceleration  of 
the  airfoil  and  the  reaction  of  the  fluid  to  this  acceleration  (Sane  2003,  Ansari  et  al.  2006). 

Wake  capturing  is  another  concept  which  is  widely  investigated  in  flapping  motion  studies.  During  reversal, 
airfoil  enters  to  the  trace  of  the  vortices  generated  during  the  previous  stroke  and  encounters  to  their  velocity 
and  acceleration  fields,  thus  resulting  in  higher  aerodynamic  forces  immediately  following  stroke  reversal 
(Sane  2003).  This  is  the  explanation  of  the  second  peak  observed  at  the  instantaneous  aerodynamic  force 
coefficients. 

This  study  aims  to  investigate  the  flow  around  a  flapping  airfoil  for  different  angles  of  attack  by  means  of 
particle  image  velocimetry  (PIV)  technique  at  Re=1000.  The  objective  is  to  show  the  influence  of  flapping 
motion  in  a  flow  where  there  exist  abundant  vortex  stmctures.  The  instantaneous  circulations  are  calculated 
for  the  vortices  generated  at  the  leading  edge  and  trailing  edge  of  the  airfoil  during  the  impulsive  start  of  the 
flapping  motion  and  during  the  pseudo-periodic  motion  where  the  influence  of  the  impulsive  start  disappears. 


2.  Description  of  Motion  and  Experimental  Setup 

The  model  is  simplified  by  use  of  a  symmetrical  airfoil  (NACA  0012).  The  flapping  motion  is  composed  by 
the  superposition  of  a  translational  motion  and  a  rotational  motion  around  a  center  of  rotation.  During  the 
translational  phase,  the  airfoil  translates  with  a  constant  velocity  until  the  time  tv  and  the  position  xv,  where  a 
rotational  motion  around  a  point  on  the  chordline  is  superposed  to  the  translational  motion  at  a  predefined 
time  ta  and  position  xa.  The  translational  velocity  is  zero  (V  =  0)  at  the  end  of  the  strokes.  When  the  time  ta 
and  the  position  is  reached,  the  airfoil  starts  to  rotate  around  its  center  of  rotation  to  have  90°  angle  of 
attack  at  the  end  of  the  stroke.  The  rotation  is  such  that  the  leading  edge  stays  as  leading  edge  during  all 
phases  of  the  motion.  The  translational  velocity  (V)  and  the  angular  velocity  (co)  variation  are  given  by  Eq.  1 
and  Eq.  2.  The  continuity  of  the  velocities  is  conserved  by  these  suggested  motions. 
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Fo  is  the  constant  translational  velocity,  a0  is  the  constant  angle  of  attack  and  T  is  the  period  of  the  motion. 

The  experimental  setup  is  a  1.5  mx  1  mx  1  m  water  tank  made  of  altuglas  and  filled  with  water  (Figure  41). 
A  rectangular  wing  of  6  cm  chord  and  50  cm  length  having  a  NACA  0012  airfoil  section  is  displaced  inside 
the  tank  by  associating  a  rotational  and  a  translational  motion.  The  center  of  rotation  is  at  the  quarter  chord 
location.  The  wing  is  delimited  by  two  rectangular  plates  made  of  epoxy  with  50  cm  x  90  cm  dimensions  in 
order  to  obtain  a  2D  flow.  PIV  measurements  are  performed  with  a  30  mJ  Quantel  Twin  Ultra  Nd:YAG  laser 
using  two  CCD  cameras  in  order  to  visualize  the  whole  flapping  domain  with  a  field  of  view  of  25  cm  x  50 
cm  with  each  one  recording  80  double  frames  per  period.  The  laser  is  coupled  with  a  spherical  lens  followed 
by  a  cylindrical  lens.  A  mirror  is  located  at  the  bottom  center  of  the  water  tank  to  spread  the  light  of  the  laser 
to  the  whole  flow  domain  as  wide  as  possible  and  to  limit  the  shadow  generated  by  the  model.  Two  CCD 
cameras  with  60  mm  objectives  and  F  number  of  4  are  used  to  record  the  entire  flow  domain. 


Figure  41  Experimental  setup  in  water  tank 


3.  Results 

The  flow  periodicity  is  investigated  for  a  single  case  chosen  with  starting  angle  of  attack  a0=60°  and 
impulsive  period  and  the  following  periods  are  compared  until  the  flow  periodicity  occurs.  The  impulsive  start 
effect  is  observed  to  disappear  before  the  7th  period.  Thereafter,  the  influence  of  different  parameters  on 
vortex  topology  and  circulation  of  the  vortices  are  represented. 
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3.1  Vortex  recognition  method 

The  velocity  gradient  is  made  of  two  parts: 
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where  i,j  (=1,2)  are  free  indices,  u;  is  corresponding  velocity  component  and  x,  is  the  corresponding 
space  coordinate  in  Cartesian  system.  The  velocity  gradient  is  summation  of  the  symmetrical  rate-of- 
strain  (deformation)  tensor  S;j  and  the  skew-symmetrical  rate-or  rotation  tensor  Q,,). 


\ 


du .  du. 


Q .  =  - 

lJ  2 


dut 

dx. 


d uj 

dx. 


(4) 


(5) 


The  second  invariant  of  velocity  gradient  (also  called  as  second  invariant  of  the  mean  rate-of- 
displacement  tensor)  is  given  by  Eq.6. 
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Nondimensional  Q  is  defined  as: 
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where  c  is  the  chord  length  and  Vo  is  the  maximum  velocity  of  the  airfoil  during  the  flapping  motion. 


For  the  2D  motion: 
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A  technique  by  use  of  Q  criteria  is  used  to  distinguish  the  instantaneous  and  local  vortex  formation.  The 
centers  of  the  instantaneous  local  vortices  are  found  by  using  maximum  local  Q  values  in  the  velocity 
flowfield.  In  addition,  positive  Q  contours  are  used  to  detect  the  boundaries  of  each  local  vortex. 


Circulation  of  each  vortex  is  calculated  by  Eq.  9. 
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where  the  integration  is  carried  out  along  the  detected  vortex  boundary.  V  denotes  the  velocity  vector  along 
the  contour  of  the  vortex  and  dl  stands  for  the  elemental  length  along  the  contour. 
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The  airfoil  starts  from  the  rest  its  motion  with  a  90°  angle  of  attack  with  zero  translational  velocity  and 
reaches  at  a  constant  angle  of  attack  a0  at  and  at  a  constant  velocity  V0  at  xv=2c.  The  rotation  and  translation 
of  the  airfoil  are  referenced  with  respect  to  the  l/4c  centre  of  rotation. 

The  flow  periodicity  is  investigated  with  starting  angle  of  attack  a0=60°  and  impulsive  period  and  the 
following  periods  are  compared  until  the  flow  periodicity  occurs.  The  impulsive  start  effect  is  observed  to 
disappear  before  the  7th  period.  Figure  42  represents  instantaneous  vorticity  distribution  beginning  from  the 
impulsive  start  of  flapping  airfoil  for  a0=60°  case.  Consecutive  PIV  images  are  represented  with  a  time  step 
of  At*=0.0125  (80  photo  per  period)  from  the  beginning  to  the  half-amplitude  location  of  the  upstroke. 

Erreur  !  Source  du  renvoi  introuvable.  represents  the  vortices  distinguished  by  the  Q  criteria  from  the 
PIV  flowfield.  The  centers  of  the  vortices  are  represented  by  the  red  symbols  for  the  counter-clockwise 
vortices  and  by  blue  symbols  for  the  clockwise  vortices. 
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Figure  42  Instantaneous  vorticity  for  a0=60°  beginning  from  impulsive  start  (photos  are  consecutive  PIV  images 
during  1st  stroke. 
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Figure  43  Instantaneous  vorticity  for  a0=60°  beginning  from  impulsive  start  (photos  are  consecutive  PIV 
images  during  1st  stroke  (The  centers  of  the  vortices  are  represented  by  the  symbols  o  for  the  counter¬ 
clockwise  vortices  and  +  for  the  clockwise  vortices). 

The  instantaneous  circulation  data  for  each  vortex  formation  during  the  impulsive  start  (1st  period)  is 
represented  in  Figure  4.  Impulsive  start  effect  is  visiple  on  the  firstly  generated  trailing  edge  vortex  TV1.  It 
has  a  negative  circulation  sign.  Both  the  generated  leading  edge  vortex  (LEVI)  and  trailing  edge  vortex  - 
TV1)  have  influences  on  the  flowfield  domain  during  the  both  upstroke  and  the  downstroke  of  the  motion. 


Figure  44  Instantaneous  non-dimensional  circulation  for  a0=60°  beginning  from  impulsive  start.  The  vortices 
generated  during  the  first  half-stroke.  (LEV:  leading  edge  vortex,  TV:  Trailing  edge  vortex,  RSV:  Rotational 
stopping  vortex). 

At  stroke  reversal,  the  shedding  of  the  leading  edge  vortex  and  rotational  vortex  forms  a  counter¬ 
rotating  pair  which  imposes  a  fluid  flow  towards  the  lower  surface  of  the  wing  at  the  onset  of  the 
next  stroke  (Birch  and  Dickinson,  2003).  The  vorticity  shedding  from  the  previous  stroke  influences 
the  force  production  (Birch  and  Dickinson  2003). 

Figure  45  represents  the  vortices  generated  during  the  first  half-stroke  of  the  flapping  motion.  As  it  is  notices, 
all  the  vortices  are  generated  before  t*=0.5.  The  LEVI  (leading  edge  vortex),  represented  with  black  circle  on 
the  figure  is  the  most  dominant  vortex  in  the  flowfield.  Until  the  end  of  the  translation  phase,  we  observed  that 
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its  circulation  increases  gradually.  When  the  airfoil  starts  to  its  rotation,  first  we  observed  that  the  circulation 
of  the  LEV  1  decreases  and  then  increases  very  strongly  until  the  end  of  the  rotational  motion. 


Figure  45  Instantaneous  vorticity  and  velocity  vectors  relative  to  flapping  translational  velocity  for  a0=60° 
(beginning  of  upstroke  during  different  periods  [from  IT  to  7T]) 


The  instantaneous  vorticity  contours  during  the  impulsive  start  are  represented  in  Figure  46  for  two  angles  of 
attack  distribution  where  the  constant  angle  of  attack  a0  is  reached  at  xa=2c. 
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3)  ccq—30°  b)  cxq~30° 


Figure  46  Instantaneous  non-dimensional  vorticity  contours  for  impulsive  start  of  the  flapping  motion  from  90° 
to  30°  (left  column)  and  from  90°  to  60°  (right  column)  by  PIV  measurements 


With  the  impulsive  rotation  and  translation  of  the  airfoil,  a  counter-clockwise  trailing  edge  vortex  is  observed 
with  a  positive  circulation.  This  trailing  edge  vortex  detaches  from  the  airfoil  surface  and  is  pushed  outside  the 
near  wake  where  takes  place  one  negative  circulation  trailing  edge  vortex.  It  is  noted  that,  this  first  impulsive 
trailing  edge  vortex  is  stronger  for  a0=30°  since  the  instantaneous  angular  velocities  are  much  bigger  for  this 
case.  In  the  mean  time,  a  leading  edge  vortex  is  growing  with  a  positive  circulation  value.  At  t*=0.0133 ,  just 
after  the  motion  of  the  airfoil  starts,  this  counter-clockwise  leading  edge  vortex  starts  to  be  visualized  clearly 
and  it  grows  up  with  the  airfoil  motion.  The  non-dimensional  circulation,  r/Voc ,  of  this  leading  edge  vortex 
for  a0=60°  case  is  nearly  twice  of  the  one  for  a0=30°  (Figure  47).  Similar  vortex  topologies  are  observed 
both  for  a0=30°  and  a0=60°  angles  of  attack  cases  however  the  traces  of  the  leading  and  trailing  edge  vortices 
are  more  attached  to  the  airfoil  surface  for  a0=30°  case. 
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Figure  47  Instantaneous  circulation  of  the  vortices  by  PIV  measurements  for  a0=30°  (left  column)  and  a0=60° 

(right  column)  at  t*=0.1 008 


Conclusion  and  Perspective 

The  aim  of  the  present  study  was  to  identify  the  unsteady  vortex  structures  during  a  flapping  motion  in  hover 
mode  at  Re=1000.  Particle  Image  Velocimetry  measurements  are  used  in  order  to  identify  vortices  generated 
during  the  flapping  motion  for  different  parameters.  The  complex  instantaneous  vortex  dynamics  is 
investigated  first  at  the  impulsive  start  region  and  then  in  the  quasi-periodic  region  where  the  impulsive  effect 
is  lost.  The  instantaneous  flow  topology  shows  the  importance  of  the  interaction  of  the  vortices  generated  with 
the  trace  of  the  vortices  from  the  previous  stroke. 
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2.2  Design  of  a  Four  Wing  Hovering  MAV  with  Cooperating  Slider 
Crank  Mechanism 

(by  Mehmet  Kilic,  M.Sc  student  ME  METU,  Asst.  Prof.  Dr.  D.  Funda  Kurtulus  AE  METU,  Asst.  Prof.  Dr. 
Yigit  Yazicioglu  ME  METU) 


This  study  examines  the  properties  of  the  two  dimensional  four  bar  mechanism  and  the  cooperating 
slider  crank  mechanisms  for  flapping  wings.  A  detailed  derivation  procedure  is  given  for  the 
cooperating  slider  crank  mechanisms.  It  is  also  noted  that  when  the  cranks  of  the  slider  crank 
mechanisms  coincide,  the  kinematics  of  the  mechanism  simplifies  greatly.  It  is  also  noted  that  with  a 
cooperating  slider  crank  mechanisms,  symmetrically  working  wings  can  easily  be  obtained.  The 
parts  of  the  micro  air  vehicle  utilizing  this  mechanism  was  designed,  manufactured  and  tested  for 
the  required  motion.  As  a  future  work,  the  flexible  membrane  wings  will  be  added  to  the  mechanism 
and  dynamic  lift  forces  will  be  measured. 


1  -  Introduction 

With  a  more  general  sense,  the  wing  of  a  micro  air  vehicle  can  make  at  most  three  rotations.  So 
that,  any  motion  required  from  the  wing  can  be  obtained.  But  an  increased  degree  of  freedom 
results  in  a  complicated  mechanism  and  the  weight  of  the  micro  air  vehicle  increases  which  is  not 
favored.  Then  researchers  turn  their  faces  to  low  degree  of  freedom  and  low  weight  mechanisms. 
To  obtain  a  one  degree  of  freedom  flapping  motion,  the  first  mechanisms  visited  are  the  two 
dimensional  and  the  three  dimensional  four  bar  mechanisms.  There  are  several  successful 
examples  on  this  line  at  the  literature  (See  Figure  48  and  Figure  49). 


Figure  48  The  transmission  designs  with  two  dimensional  four  bar  mechanisms. 
Microbat  at  the  left  [1,  2],  and  Mechanical  Humming  Bird  at  the  right  [3] 
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Figure  49  The  transmission  designs  with  two  dimensional  four  bar  mechanisms.  Delfly  I  [4]  at  the  up 
The  transmission  designs  with  three  dimensional  four  bar  mechanisms.  Delfly  II  [4]  at  the  down 

2  -  Design  of  a  Two  Dimensional  Four  Bar  Mechanism  for  Flapping 

With  a  four-bar  mechanism,  more  specifically  with  a  crank  rocker  type  one,  we  can  easily  obtain  a 
flap  motion  (See  Figure  50). 


Figure  50  Crank  rocker  type  four  bar  mechanism.  AB  is  the  crank  and  DC  is  the  rocker. 


At  hovering  flight,  we  expect  the  wing  make  the  same  motion  during  the  downstroke  and  the 
upstroke.  We  cannot  totally  fulfill  this  task  with  a  four  bar  mechanism,  but  the  timing  of  the 
downstroke  and  the  upstroke  can  be  aligned.  When  9  is  0  and  180°,  the  rocker  can  happen  to  be  at 
the  dead  centers.  These  two  positions  are  shown  at  Figure  50.  If  the  rocker  is  forced  to  swing  Ay/ 
much  during  the  flaps,  then  there  appear  the  following  relations. 


\DC\ 


M 

sin(A^/2) 


\AD\ 


\ab\ 


tan(A^/2) 


(1) 

(2) 


62 


Notice  here  that  we  can  play  with  three  parameters:  \AB\ ,  \BC\  and  A y/.  For  \AB\  =  10  mm,  \BC\  = 

50  mm  and  Ay/  =  80  degrees,  the  deviation  of  y/  with  respect  to  6  is  given  in  Figure  55  (.  For 
comparison,  a  cosine  curve  is  also  given.  As  it  is  seen,  the  resultant  swinging  motion  deviates  from 
cosine  curve.  When  the  crank  is  rotating  with  a  constant  angular  velocity,  the  downstroke  and  the 

upstroke  takes  same  amount  of  time  but  the  motion  differs  from  each  other  greatly.  If  \BC\  were 

selected  to  be  longer,  the  motion  curve  would  approach  the  cosine  curve  more. 

3  -  Design  of  a  Cooperating  Slider  Crank  Mechanism  for  Flapping 

Consider  two  inline  slider  crank  mechanisms  using  the  same  slider  (See  Figure  51).  The  link  AB 
drives  the  blue  slider  crank  mechanism.  Then,  the  slider  drives  the  pink  slider  crank  mechanism. 
Although  link  AB  makes  full  rotation,  now  link  ED  makes  a  swinging  motion. 


We  want  to  correlate  the  angular  position  of  link  ED  with  the  angular  position  of  link  AB.  That  is,  a 
correlation  between  a  and  6 .  An  arbitrary  function  relating  a  to  6  cannot  be  realised  for  all 
angles  to  be  satisfied  with  the  cooperating  slider  crank  mechanisms,  yet,  we  can  approximate  it 
satisfying  the  function  at  the  precision  points.  Like,  when  0  is  9X ,  a  is  ax  and  when  0  is  02 ,  a  is 

a2 ,  etc. 

When  we  assign  the  dimensions  of  links  AB  and  BC,  we  will  know  the  locations  of  the  slider  point  C. 
Let  the  distance  between  the  slider  and  point  E,  (that  is,  \EC\ ),  is  S.  Now  the  correlation  problem 
reduces  to  correlating  the  rotations  of  the  link  ED  with  the  slider  distance  S. 


The  possible  three  positions  of  the  slider  crank  mechanism  are  given  in  Figure  52.  At  the  first 
position,  to  fully  define  the  mechanism,  we  have  in  total  4  parameters.  Namely,  \ED\ ,  ax ,  \DC\  and 
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j3x.  On  the  other  hand,  satisfying  the  first  position  imposes  only  two  kinematic  equations.  That  is,  4 
parameters  and  corresponding  2  equations.  Now  we  assume  the  task  as  follows;  when  the  slider 
translates  from  Sx  to  S2 ,  link  ED  deviates  a2  much.  Here,  we  gain  another  two  kinematic  equations 

with  the  cost  of  a  new  parameter,  J32 .  Then  in  total,  5  parameters  and  corresponding  4  equations.  In 
a  similar  manner,  satisfying  a  third  position  similar  to  the  second,  we  will  have  6  parameters;  \ED\ , 

ax ,  \DC\,  Px,  p2  and  J33,  and  the  corresponding  6  kinematic  equations.  This  problem  can  be 
solved  analytically  and  it  is  explained  below. 

In  complex  plane,  the  6  kinematic  equations  can  be  expressed  with  3  vector  equations, 


where 


ED  +  DC  =  Sx 

(3) 

ED  e'“2  +DCeiP2  =  S2 

(4) 

~ED  e'“3  +  ~DC  ei/h  =  S3 

(5) 

ED  =  \ED\  eia' 

(6) 

DC  =  \DC\eip' 

(7) 

Observing  the  equations  3,  4  and  5;  we  see  that  they  are  linear  in  terms  of  ED  and  DC  ,  but  the 
equation  set  is  over  defined.  Then,  the  determinant  of  the  augmented  matrix  must  be  zero. 


Manipulating  the  determinant, 


where 


-A,  +  A2  eiPl  -  A3  eiPl  =  0 

A \  =  S3eia>  -S2eia> 

A ~  =  S3-Sxeia> 

A ]  =  S2-Sxeia> 


(8) 


(9) 

(10) 
(11) 
(12) 


The  solution  of  equation  9  can  be  treated  graphically  (See  Figure  53).  It  is  clear  that  there  are  either 
two  solution  set  or  no  solution  (remember  that  the  dashed  circles  can  not  intersect). 
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Figure  53  Graphical  resolution  to  Equation  9 


After  finding  the  value  of  /?2 ,  we  insert  it  into  Equation  4,  and  Equations  3  and  4  can  now  be  solved 

for  ED  and  DC  .  Then  using  Equations  6  and  7,  we  calculate  the  remaining  unknown  parameters 
and  the  mechanism  is  fully  designed. 

With  the  four  bar  mechanism  that  have  been  designed,  we  guaranteed  that,  y/  is  40°  when  6  is  0 
and  y/  is  -40  0  when  6  is  180  °.  With  a  cooperating  slider  crank  mechanisms,  we  can  further  force 
the  system  to  y/  is  0  when  0  is  90  °.  So  the  y/  deviation  with  respect  to  6  will  imitate  the  cosine 
curve  better.  The  problem  was  solved  for  |^45|=10,  \BC\  =40,  =15,  02=  90°,  a2  =  40°,  02  -  180° 

and  a3  =  80°.  One  of  the  two  possible  solution  sets  resulted  in  infinite  \ED\  and  \DC\ ,  then  only  one 
solution  was  possible  and  shown  in  Figure  54. 


i4fl|=10,  \BC\  =40,  |^|=15,  \ED\  =15,94  and  \DC\  =57,43 


Cooperating  slider  crank  mechanism  is  superior  to  two  dimensional  four  bar  mechanism  in  some 
other  ways  other  than  approximating  the  cosine  curve  better.  Firstly,  it  benefits  from  the  symmetry  of 
the  inline  nature.  By  default,  the  ^deviation  with  respect  to  6  will  be  symmetric  with  respect  to  6  = 
180°  line  like  the  cosine  curve.  That  is,  the  upstroke  of  the  flap  is  totally  the  same  as  the 
downstroke.  The  two  dimensional  four  bar  mechanism,  on  the  other  hand,  does  not  have  such  a 
property  (See  Figure  55).  Secondly,  one  can  easily  be  place  a  second  flapping  wing  to  the 
cooperating  slider  crank  mechanism  which  flaps  symmetrical  to  the  first  wing  (See  Figure  56). 
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Figure  56  Cooperating  slider  crank  mechanisms  actuating  2  wings. 

ED  and  ED’  are  symmetrically  working  wing  arms. 

4  -  The  Cooperating  Slider  Crank  Mechanism  used  in  our  design 

After  feeling  the  merits  of  the  cooperating  slider  crank  mechanism,  we  intended  to  use  it  for  the 
proposed  design.  There  were  some  restrictions  for  the  parts.  First  of  all,  we  desire  to  obtain  the  wing 
wing  interaction  at  the  ends  of  the  two  strokes.  Therefore,  the  only  way  to  rotate  the  crank  is  driving 

it  with  a  shaft  passing  through  a  hole  at  the  rotation  center  of  the  wings.  That  means  the  \AE\  will  be 

zero.  This  special  case  resulted  in  a  more  facilitated  design  procedure  for  the  cooperating  slider 
crank  mechanism.  Instead  of  the  design  procedure  given  in  part  3,  the  following  rules  suffice. 

When  \AD\,  |Z)C| ,  Ay/  are  given, 


AB\  =  \AD\  sin(A^/2) 


(13) 
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The  problem  was  solved  for  |^4Z>|  =25,  \DC\  =45,25  and  A y/  =80°  and  shown  in  Figure  57. 


Figure  57  Cooperating  slider  crank  mechanisms  satisfying  the  precision  points. 
\AD\  =25,  \DC\  =45,25,  Ay/  =80°  \  AB\  =16,07and  \BC\  =41 


The  Equation  14  guaranties  that  line  DB  is  always  perpendicular  to  line  AC.  Then  a  simple  explicit 
relation  between  y/  and  6  can  be  obtained. 


y/  =  a  sin 


r\AB\ 


cos(#) 

J 


(15) 


The  mechanism  working  with  this  principle  was  manufactured  and  tested  for  the  required  motion 
(See  Figure  58).  All  parts  were  made  up  of  aluminum.  The  wing  spars  are  2  mm  thick  carbon  fiber 
rods. 


Currently,  we  work  on  the  membrane  wing  geometry  and  the  wing  stiffener  locations.  We  intend  to 
measure  real  time  lifting  force  history  using  sensitive  force  transducers.  Hopefully,  the  correct 
dimensions  will  be  found  for  the  most  efficient  lifting  force  generation. 
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Figure  58  The  mechanical  flapper  using  cooperating  slider  crank  mechanisms. 
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2.3  Experimental  Setups  at  AE  METU 


(by  Hasan  Hizli,  M.Sc  student  AE  METU,  Asst.  Prof.  Dr.  D.  Funda  Kurtulus  AE  METU  with  collaboration  of  Robotics 
&  Mechatronics  Technologies  Ltd) 


As  it  is  explained  in  Report  1,  we  were  at  the  stage  of  establishing  the  new  experimental  setup  and  a 
water  supply  system  to  the  laboratory  that  is  given  to  me  as  a  research  area.  The  water  supply  system 
has  been  added  and  the  new  experimental  setup  has  been  designed  for  water  tank  (Figure  59).  Tests  of 
the  experimental  setup  will  be  started  on  September  2009. 


Figure  59  Water  tank  (2m  length) 


Introduction 

Research  Pool  Positioning  System  (RPPS)  is  an  advanced  laboratory  platform  designed  to  accurately 
position  airfoils  within  research  pools  during  experiments.  The  system  enables  users  to  move  and 
locate  airfoil  equipment  in  two  dimensions  both  longitudinal  and  rotational  and  provides  logging  for 
the  complete  motion  during  operation.  RPPS  can  control  and  log  the  motion  of  the  equipment  in  real¬ 
time  through  a  PC  installed  Data  Acquisition  Board  and  Matlab/Simulink™  software. 

RPPS  can  perform  the  following  demonstrations  and  experiments: 

Linear  position  and  rotation  control  of  the  cart  mounted  airfoil. 

Real-time  linear  and  rotational  position  data  tracing  and  logging. 

On  the  fly  parameter  updating. 
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Besides  the  basic  functionalities  given  above,  the  system  can  be  equipped  with  additional  sensors  to 
make  calibration,  estimation,  identification  and  control  practices  in  the  laboratory. 


The  block  diagram  of  RPPS  is  presented  in  Figure  60. 


Control 

Computer/ 

Software 


Data 

Acquisition 

Board 


Power  Supply 


Linear  and  Rotational  Positioning  Unit 


Figure  60  System  Block  Diagram. 


System  Components 

System  shown  in  Figure  60  is  explained  in  details  below. 

Positioning  Unit 

Positioning  Unit  (PU),  is  the  basic  mechanical  structure  which  holds  the  cart  via  the  rack  and  the 
linear  bearing.  The  main  construction  is  composed  of  high  quality  machined  aluminum  housing  and 
holds  a  custom  precision  machined  rack.  A  picture  of  PU  is  available  in  Figure  61. 
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Figure  61  Positioning  unit. 


Cart 

The  Cart  carries  the  airfoil  equipment  and  moves  along  the  rack  via  the  pinion  mounted  on  the 
motor.  The  linear  position  control  of  the  airfoil  is  achieved  with  the  motion  of  the  cart.  A  cart 
mounted  independent  mechanism  controls  the  rotational  motion  of  the  airfoil. 

The  cart  consists  of  the  carrier  platform  and  two  drive  motors  with  embedded  gear  trains  and 
encoders  for  translational  and  rotational  motion  of  the  airfoil. 


The  Cart  uses  high  quality  gear  motors  for  translational  and  rotational  motion  along  the  rack  (Figure 
62).  The  geared  motor  has  an  embedded  encoder  for  measuring  direct  position  of  the  pinion  during 
cart  motion.  The  precision  gear  motor  is  a  low  power  motor  with  high  efficiency  for  smooth  and 
quite  operation  in  a  small  volume.  The  specifications  for  the  motor  are  given  below. 


Figure  62  Cart  motor  with  embedded  encoder. 
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Motor: 


Motor  Data 

Units  Value 

Assigned  power  rating 

W  5.0 

Nominal  voltage 

Volt  6.0 

Gearhead: 


Gearhead  Data 

Units  Value 

Planetary  gearhead  reduction 

19:1 

Number  of  stages 

2 

Digital  Encoder: 

Magnetic  digital  encoder  provides  16  counts  per  turn.  The  encoder  is  embedded  to  the  drive  motor. 

Interface  Box 

Interface  Box  (IB)  is  the  main  junction  box  between  the  data  acquisition  board,  power  supply  and 
motors.  The  block  diagram  of  the  box  is  given  below: 
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Figure  63  Interface  Box  block  diagram. 


Control  inputs  to  the  motor  and  sensor  signals  from  the  encoder  are  passed  through  the  interface  box. 
IB  also  contains  electronics  that  provides  motor  driving  and  power  distribution  among  units. 

Motor  Driver  Electronics 

Positioning  unit  driving  motors  is  driven  by  pulse  width  modulation  (PWM)  signals  through 
commands  received  from  the  controller  PC  through  data  acquisition  board  (Figure  64). 


Figure  64  Motor  driver  electronics  (Robotsan  RS-MD02  and  rsPIC  units). 
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Motor  driving  is  accomplished  through  Robotsan’s  RS-MD02  Motor  Driver  and  rsPIC  Intelligent 
Controller  units. 

Data  Acquisition  Board 

Data  Acquisition  Board  (DAQ  Board)  serves  as  the  signal  acquisition  and  command  generation 
electronics  onboard  the  PC.  Encoder  signals  and  motor  driving  logic  and  commands  are  generated  on 
the  DAQ  Board.  The  board  is  also  compatible  with  Matlab/Simulink™  and  Real-Time  Workshop® 
software  enabling  real-time  implementations. 

Data  Acquisition  Board  used  within  the  system  is  shown  in  Figure  65. 


Figure  65  Data  Acquisition  Board. 
Features  of  the  board  are  listed  below: 

•  Sixteen  16-bit  differential  A/D  inputs,  15  kHz  rate 

•  Four  14-bit  D/A  outputs,  20  kHz  update  rate 

•  48  bi-directional  digital  signals 

•  Six  quadrature  encoder  inputs 

•  Driver  for  Real-Time  Windows  Target 

•  Driver  for  xPC  Target 

•  Driver  for  Windows 


Real  Time  Control  Software 

Real-Time  Windows  Target™  rapid  prototyping  software  is  a  PC  solution  for  prototyping  and 
testing  real-time  systems.  Real-Time  Windows  Target  software  uses  a  single  computer  as  a  host  and 
target.  On  this  computer,  you  use  the  MATLAB®  environment,  Simulink®  software,  to  create  models 
using  Simulink®  blocks  and  Stateflow®  diagrams. 
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After  creating  a  model  and  simulating  it  using  Simulink  software  in  normal  mode,  you  can  generate 
executable  code  with  Real-Time  Workshop®  code  generation  software  and  the  C/C++  compiler. 
Then  you  can  run  your  application  in  real  time  with  Simulink  external  mode. 

Integration  between  Simulink  external  mode  and  Real-Time  Windows  Target  software  allows  you  to 
use  your  Simulink  model  as  a  graphical  user  interface  for 

Signal  visualization  —  Use  the  same  Simulink  Scope  blocks  that  you  use  to  visualize  signals 
during  a  non-real-time  simulation  to  visualize  signals  while  running  a  real-time  application. 

Parameter  tuning  —  Use  the  Block  Parameter  dialog  boxes  to  change  parameters  in  your 
application  while  it  is  running  in  real  time. 

Operating  RPPS 

Simulink  Interface 

This  section  describes  the  Simulink  blocks  used  to  interface  with  the  RPPS  hardware.  The  following 
explains  these  interfaces  and  provides  notes  on  their  usage: 

Actuator  Block 

This  block  enables  the  user  to  provide  a  control  signal  to  drive  the  motors  for  both  translational  and 
rotational  motion. 

The  input  of  the  block  is  a  voltage  between  0-5V’s. 

The  outputs  of  the  blocks  are: 

-  an  analog  signal  on  the  DAC  output  of  DAQ  providing  the  required  drive  voltage  on  the  motor, 

-  control  of  two  Digital  Outputs  on  the  DAQ  providing  direction  control  of  the  motor. 


Figure  66  Actuator  Block  used  for  driving  the  motors. 
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The  block  determines  the  sign  of  the  control  signal  which  then  calculates  the  required  PWM  value 
and  the  direction  of  the  actuator.  The  user  is  only  required  to  provide  the  control  input. 

The  subsystem  view  of  the  Actuator  Block  is  presented  in  the  next  figure. 


Figure  67  Actuator  Block  details. 


Actuator  Feedback  Encoder  Block 

This  block  enables  the  user  to  read  the  incremental  position  of  the  drive  motors. 

The  output  of  this  block  is  the  motor  revolutions  in  encoder  steps.  The  Encoder  Input  block  is  used  to 
receive  the  encoder  signals  from  the  DAQ  (Channel  3).  A  view  of  the  block  and  conversion  factors 
are  given  below: 


Figure  68  Feedback  Encoder  Blocks. 
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The  block  parameters  are  also  presented  in  the  next  figure.  In  order  to  increase  the  resolution  of 
feedback  data  available  the  encoder  is  read  in  quadrature  mode  where  1024x4=4096  counts  per  turn 
is  achieved  during  encoder  readings.  The  sample  time  Ts  (i.e.,  500Hz)  should  be  set  from  within  the 
Matlab  window  prior  to  operation. 


Figure  69  Actuator  Feedback  Encoder  Input  Block  Parameters. 
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